#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"

#define ITER            25
#define PRECISION       100000000
#define BASE            10
#define DIGITS          30000000

int
main(int argc, char *argv[])
{
    int i, iter;
    fprintf(stderr,"Calculando Borwein Sequencial com %d iteracoes...\n",iter);

    /* main variables */
    mpf_t a, b, p;

    /* aux */
    mpf_t sqrt_a, inv_sqrt, sum1, sum2, res1, res2, mult, ab;

    /* Inicialization of the variables */
    mpf_set_default_prec(PRECISION);
    mpf_inits(ab, a, b, p, sqrt_a, inv_sqrt, sum1, sum2, res1, res2, mult, NULL);

    /* a = sqrt(2) */
    mpf_init_set_ui(a, 2);
    mpf_sqrt(a, a);

    /* b = 0 */
    mpf_init_set_ui(b, 0);

    /* p = 2 + sqrt(2) */
    mpf_init_set_ui(p, 2);
    mpf_add(p, p, a);

    if(argc == 2)
        iter = atoi(argv[1]);
    else iter = ITER;

    for(i = 0; i < iter; i++)
    {
        mpf_set(ab, a);

        /* an1 = (sqrt(an) + (1 / sqrt(an))) / 2; */
        mpf_sqrt(sqrt_a, a);
        mpf_ui_div(inv_sqrt, 1, sqrt_a);
        mpf_add(sum1, sqrt_a, inv_sqrt);
        mpf_div_ui(a, sum1, 2);

        /* bn1 = ((1 + bn) * sqrt(an)) / (an + bn) */
        mpf_add_ui(sum1, b, 1);
        mpf_sqrt(sqrt_a, ab);
        mpf_mul(mult, sum1, sqrt_a);
        mpf_add(sum2, ab, b);
        mpf_div(b, mult, sum2);

        /* p = ((1 + a) * p * b) / (1 + b) */
        mpf_add_ui(sum1, a, 1);
        mpf_mul(res1, sum1, p);
        mpf_mul(res2, res1, b);
        mpf_add_ui(sum2, b, 1);
        mpf_div(p, res2, sum2);
    }

    /* Cleaning the memory */
    mpf_out_str(stdout, BASE, DIGITS, p);
    mpf_clears(ab, a, b, p, sqrt_a, inv_sqrt, sum1, sum2, res1, res2, mult, NULL);

    return 0;
}
